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Abstract 

In this paper, a non-uniform rational B-spline based iso-geometric finite element method is used to study the 
static and dynamic characteristics of functionally graded material (FGM) plates. The material properties are 
assumed to be graded only in the thickness direction and the effective properties are computed either using 
the rule of mixtures or by Mori-Tanaka homogenization scheme. The plate kinematics is based on the first 
order shear deformation plate theory (FSDT). The shear correction factors are evaluated employing the energy 



^ ■ equivalence principle and a simple modification to the shear correction factor is presented to alleviate shear 

o ■ 

: locking. Static bending, mechanical and thermal buckling, linear free flexural vibration and supersonic flutter 
analysis of FGM plates are numerically studied. The accuracy of the present formulation is validated against 



available three-dimensional solutions. A detailed numerical study is carried out to examine the influence of 
the gradient index, the plate aspect ratio and the plate thickness on the global response of functionally graded 
material plates. 
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1. Introduction 

Since its introduction to decrease the thermal stresses in propulsion systems and in airframes for space appli- 
cation functionally graded materials (FGMs) have led researchers to investigate the structural behaviour of 
such structures. FGMs are considered to be an alternative for certain class of aerospace structures exposed to 
high temperature environment. FGMs are characterized by a smooth transition from one material to another, 
thus circumventing high inter-laminar shear stresses and de-lamination that persists in laminated composites. 
Thus, for structural integrity, FGMs have advantages over the fiber-matrix composites. 

1.1. Background 

The investigation of the static and the dynamic behaviour of FGM structures is fairly well covered in the 



ana. 



literature. Some of the important contributions are discussed here. Different plate theories, viz, FSDT 
second and other higher order accurate theory have been used to describe plate kinematics. Existing 

approaches in the literature to study plate and shell structures made up of FGMs uses finite element method 
based on Lagrange basis functions 2, 8, 4], meshfree methods [f], 0]. All existing approaches show shear locking 
when applied to thin plates. Different techniques by which the locking phenomenon can be suppressed include: 

• Retain the original interpolations and subsequently use an optimal integration rule for evaluating the 
bending and the shear terms; 



Mixed interpolation technique [9(]; 

Use field redistributed substitute shape functions 



• Discrete shear gap method [121 ]; 
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• Stabilized conforming nodal integration 



H i.e., 



strain smoothing, SFEM 



Use p-adaptivity, for example Moving Least Square approximations 

n 



a. 



He et al. , [17J presented a finite element formulation based on thin plate theory for the vibration control of FGM 



plates with integrated piezoelectric sensors and actuators under mechanical load, whereas Liew et al., 



181 ] have 



analyzed the active vibration control of plates subjected to a thermal gradient using shear deformation theory. 



The parametric resonance of FGM plates is discussed in [19| by Ng et al., based on Hamilton's principle and 
the assumed mode technique. Yang and Shen j^ol . have analyzed the dynamic response of thin FGM plates 
subjected to impulsive loads using a Galerkin Procedure coupled with modal superposition methods, whereas, 
by neglecting the heat conduction effect, such plates and panels in thermal environments have been examined 
based on shear deformation with temperature dependent material properties |3Jl. The static deformation and 

n 

vibration of FGM plates based on higher-order shear deformation theory is studied by Qian et al., [5J using the 
meshless local Petrov- Galerkin method (MLPG) and Natarajan and Ganapathi 3] using shear flexible elements. 

2l| presented analytical solutions for simply supported rectangular FGM plates based on second- 



Matsunaga 



22 



231 ] for vibrations of 



order shear deformation theory, whereas, three dimensional solutions are proposed in 
simply supported rectangular FGM plates. Reddy presented finite element solutions fc 

of FGM plates and Ferreira et al., [g] performed dynamic analysis of FGM plates based on higher order shear and 

24f | and Javaheri and Eslami 



normal deformable plate theory using MLPG. Birman 



251 ] have studied buckling 



261 ] analyzed the thermo-mechanical 



of FGM plates subjected to in-plane compressive loading. Woo et al., 
postbuckling behaviour of plates and shallow cylindrical FGM panels using a classical theory. Ganapathi et 
al., 8(], using a C° shear flexible quadrilateral element, studied buckling of non-rectangular FGM plates under 



mechanical and thermal loads. Prakash and Ganapathi 
panels exposed to supersonic flow. Haddadpour et al 



271 ] studied the linear f 



281 ] and Sohn and Kim 



29 



utter characteristics of FGM 



301 ] investigated the nonlinear 



aspects of flutter characteristics using the finite element method. FGM 



develop flaws. Recently, Yang and Chen 



311 ] and Kitipornchai et al., 



ates, like other plate structures, may 



321 ] studied the dynamic characteristics 



of FGM beams with an edge crack. Natarajan et al, 33j, [34| and Baiz et al., 35( studied the influence of the 
crack length on the free flexural vibrations of FGM plates using the XFEM and smoothed XFEM, respectively. 



1.2. Approach 

The main objective of this paper is to investigate the potential of NURBS based iso-geometric finite element 
methods to study the static and dynamic characteristics of Reissner-Mindlin plates. The present formulation 
also suffers from shear locking when lower order NURBS functions are used as basis functions, da Vaiga et 
al., 36(] showed that the shear locking phenomena can be suppressed by using higher order NURBS functions. 
A similar approach was employed to suppress shear locking in the element-free Galerkin method 3^1 . In this 
paper, we propose a simple technique to suppress shear locking, which relies on the introduction of an artificial 



shear correction factor 



381 ] when lower order NURBS basis functions are used. The drawback of this approach 



is that the shear correction factor is problem dependent. 
1.3. Outline 

The paper is organized as follows. A brief overview on functionally graded materials and Reissner-Mindlin plate 
theory is presented in the next section. Section [3] presents an overview of NURBS basis functions and a simple 
correction to the shear terms to alleviate shear locking. The efficiency of the present formulation, numerical 
results and parametric studies are presented in Section followed by concluding remarks in the last section. 



2. Theoretical Formulation 

2.1. Functionally graded material 

A rectangular plate made of a mixture of ceramic and metal is considered with the coordinates x, y along the in- 
plane directions and z along the thickness direction (see Figure ([T|)). The material on the top surface (z = h/2) 
of the plate is ceramic rich and is graded to metal at the bottom surface of the plate (z = —h/2) by a power 



4 



law distribution. The effective properties of the FGM plate can be computed by using the rule of mixtures or 
by employing the Mori-Tanaka homogenization scheme. 

Let Vi(i = c, m) be the volume fraction of the phase material. The subscripts c and m refer to ceramic and 
metal phases, respectively. The volume fraction of ceramic and metal phases are related by V c + V m = 1 and V c 
is expressed as: 

where n is the volume fraction exponent (n > 0), also known as the gradient index. The variation of the 
composition of ceramic and metal is linear for n =1, the value of n = represents a fully ceramic plate and any 
other value of n yields a composite material with a smooth transition from ceramic to metal. 

Rule of mixtures 

Based on the rule of mixtures, the effective property of a FGM is computed using the following expression: 

P = P C V C + P m V m (2) 

Mori-Tanaka homogenization method 

Based on the Mori-Tanaka homogenization method, the effective Young's modulus and Poisson's ratio are 
computed from the effective bulk modulus K and the effective shear modulus G as [4] 



K C K — K m 
K c -K rn 

G c s — G m 
G c — G m 



Vc 

1 "T y m 3K m +4G m 
Vc 

(Gc — Gm) 
■ (Gm + fl) 



1 + K 



(3) 
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where 



, _ G m (9K m + 8G 

11 ~ 6(K m + 2G m ) [ ) 

The effective Young's modulus E e g and Poisson's ratio v e s can be computed from the following relations: 



9K c gG c K 

-C'eff 



3KeK + G e fi 



UcS ~ 2(3K cS + G cff ) (5) 

The effective mass density p is computed using the rule of mixtures. The effective heat conductivity K e g and 
the coefficient of thermal expansion a e g is given by: 



K cS Kr. 



Vr 



K 



1 + v m 



(Kg ^m) 
3 K-/j). 



a eff - a r 



l 



a, 



\K C K m ) 



(6) 



Temperature distribution through the thickness 

The temperature variation is assumed to occur in the thickness direction only and the temperature field is 
considered to be constant in the xy-plane. In such a case, the temperature distribution along the thickness can 
be obtained by solving a steady state heat transfer problem: 



d 

dz 



dT 



0, T = T c at z = h/2; T = T m at z = -h/2 



(7) 
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The solution of Equation ([7]) is obtained by means of a polynomial series [39( as 



T(z) = T m + (T c - T m )7]( Zl h) 



(8) 



where, 



V(z,h) = — 



2z + h 
2h 



(2n + 



2z + h 



{An + 1)«4 ^ 2h 



2h 



2z + h 



2n+l 



(n + 1)k Ti 



(3n + 1)k 3 % 



2z + h 
2h 

2z + h 



n+l 



4n+l 



(5n + 1)k| 



2h 



2z + h 
2h 



+ 



3n+l 



5n+l 



0) 



C = 1 



+ 



(n + V)K m (2n + !)«£, (3n + 1)4 



+ 



(4n + 1)«4 (5n + 1)k 



(10) 



Reissner-Mindlin Plates 

The displacements ti, io at a point (x, y, z) in the plate (see Figure ([I])) from the medium surface are expressed 
as functions of the mid-plane displacements u ,v ,w and independent rotations 6 x ,0y of the normal in yz and 
xz planes, respectively, as: 



u(x,y,z,t) = u (x,y,t) + z9 x (x,y,t) 
v(x,y,z,t) = v (x,y,t) + z6 y (x,y,t) 

w(x,y,z,t) = w (x,y,t) (11) 
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where t is the time. 




(a) (b) 
Figure 1: (a) coordinate system of a rectangular FGM plate, (b) Coordinate system of a skew plate 

The strains in terms of mid-plane deformation can be written as: 



£ = < 



The midplane strains e p , bending strain si,, shear strain e s in Equation (|12p are written as: 









Bp 




ze b | 


H 




l°J 




, £s J 



(12) 



£« = < 



U, 



o.x 



J o,y 



£b = < 



£.s = < 




(13) 



where the subscript 'comma' represents the partial derivative with respect to the spatial coordinate succeeding 



it. The membrane stress resultants N and the bending stress resultants M can be related to the membrane 
strains, e p and bending strains Sb through the following constitutive relations: 



N 



N xx 



Nyy 

N x 



} = Ae v + Be b 



M 



mi 



> = Be p + D b e b 



M 



j-y 



(14) 



where the matrices A = Aij,~B = Bij and D b = Dij;(i,j = 1,2,6) are the extensional, bending-extensional 
coupling and bending stiffness coefficients and are defined as: 



rh/2 _ 

{Aij, B^, Dij} = / Qij{l, z, z 2 } dz 

J-h/2 



(15) 



Similarly, the transverse shear force Q = {Q X z,Qyz} is related to the transverse shear strains e s through the 
following equation: 



Qij — ^ij^ 



(16) 



where Eij = j_i^ 2 Q v i v j dz; = 4, 5) is the transverse shear stiffness coefficient, Vi, Vj is the transverse shear 
coefficient for non-uniform shear strain distribution through the plate thickness. The stiffness coefficients Q,^ 
are defined as: 
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o -o - E{z) ■ o - uE{z) • o -o -o 

Vll — ^22 — 2 _ j,2 ' ^12 ~~ I _ j,2 ' ^16 ~~ ^26 _ u 

Q44 = Q55 = ^66 = 2(T+^) 

where the modulus of elasticity E{z) and Poisson's ratio v are given by Equation ([5]). The thermal stress 
resultant N th and the moment resultant M th are: 



N th 

x ' XX 


h/2 


1 




( /V th 

yy 


> = J Q i3 a{z,T) < 


1 


> AT(z) dz 


N th 


-h/2 









h/2 


' > 

1 




yy 


> = J Qija(z,T) < 


1 


> AT(z) z dz 


M th 

A,J -xy 


-h/2 




1. > 





(18) 

where the thermal coefficient of expansion a(z, T) is given by Equation © and AT(z) = T(z) — T Q is the 
temperature rise from the reference temperature and T Q is the temperature at which there are no thermal 
strains. The strain energy function U is given by: 



U(S) = l [ {e^Ae p + e^Be b + ejBe p + ejT>e b + e^Be 8 } dtt 
1 Jn 



(19) 



where S = {u,v,w,6 x ,9 y } is the vector of the degree of freedom associated to the displacement field in a finite 
element discretization. Following the procedure given in 4(3], the strain energy function U given in Equation 
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(fT9|) can be rewritten as: 



U(S) = ^S T KS (20) 
where K is the linear stiffness matrix. The kinetic energy of the plate is given by: 

T(S) = \j [p(u 2 + vl + wl) + i(el + el)) dQ (21) 

where p = P( z ) ^z, I = f_fjj 2 z2 p( z ) dz and p(z) is the mass density that varies through the thickness of 
the plate. When the plate is subjected to a temperature field, this in turn results in in-plane stress resultants, 
N th . The external work due to the in-plane stress resultants developed in the plate under a thermal load is 
given by: 



V(S) 



+ 



24 L 



N xx {Qx,x + 8y,x) + N yy (@l,y + ^y,y) + 2N xy {&x,xOx,y + ®y,xQy,\ 



(22) 



The governing equations of motion are obtained by writing the Lagrange equations of motion given by: 



d 


~d(T-U)~ 




~d{T-U)~ 


di 


[ ddi _ 




[ ddi \ 



0, i = 1, 2, ■ ■ ■ , n 



(23) 



The governing equations obtained using the minimization of total potential energy are solved using Galerkin 
finite element method. The finite element equations thus derived are: 
Static bending: 

(K + K G ) S = F (24) 
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Free vibration: 



MS + (K + K G ) 6 = 



(25) 



Buckling analysis: 



Mechanical Buckling. 



(K + XK G )6 = 



(26) 



Thermal Buckling. 



(K + ATK G ) 6 = 



(27) 



where 6 is the vector of degree of freedom associated to the displacement field in a finite element discretization, 
AT(= T c — T m ) is the critical temperature difference, A is the critical buckling load and K, Kg are the linear 
stiffness and geometric stiffness matrices, respectively. The critical temperature difference is computed using a 
standard eigenvalue algorithm. 

Flutter analysis: The work done by the applied non-conservative loads is: 



Jn 

where Ap is the aerodynamic pressure. The aerodynamic pressure based on first-order, high Mach number 
approximation to linear potential flow is given by: 



where p a ,U a , and 9' are the free stream air density, velocity of air, Mach number and flow angle, respectively. 




(28) 




(29) 



12 



The static aerodynamic approximation for Mach numbers between v2 and 2 is 41 ] 



Ap 



PaUl 



dw , dw . / 
— — cos + — — sin 

ox oy 



(30) 



Substituting Equation (|2(jp - (|28|) in Lagrange's equations of motion, the following governing equation is obtained: 



MS + (K + XA)S = 



(31) 

— uj 2 8, the following algebraic equation is 
obtained: 



[(K + A A) -w 2 M] 6 = 



(32) 



where K is the stiffness matrix, M is the consistent mass matrix, A 



A is the aerodynamic force matrix 



and uj is the natural frequency. When A =0, the eigenvalue of ui is real and positive, since the stiffness matrix 
and mass matrix are symmetric and positive definite. However, the aerodynamic matrix A is unsymmetric and 
hence complex eigenvalues u> are expected for A > 0. As A increases monotonically from zero, two of these 
eigenvalues will approach each other and become complex conjugates. In this study, A cr is considered to be the 
value of A at which the first coalescence occurs. 



3. Non-Uniform Rational B-Splines 

In this study, the finite element approximation uses NURBS basis function. We give here only a brief introduction 

441 ] . The key ingredients in the construction 



to NURBS. More details on their use in FEM are given in (43l . 
of NURBS basis functions are: the knot vector (a non decreasing sequence of parameter values, < = 
0, 1, • • • ,m — 1), the control points, Pi, the degree of the curve p and the weight associated to a control point, 



w. The i th B-spline basis function of degree p, denoted by Ni p is defined as: 
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1 if & < £ < 6+1 
else 



&+p+l - Ci+l 



(33) 



A degree NURBS curve is defined as follows: 



c(0 



E ^, P (O^Pi 

i=0 

m 

E AUs>* 

i=0 



(34) 




Figure 2: non-uniform rational B-splines, order of the curve = 3 
where Pj are the control points and Wi are the associated weights. Figure (|2|) shows the third order non-uniform 
rational B-splines for a knot vector, E = {0, 0, 0, 0, 1/3, 1/3, 1/3, 1/2, 2/3, 1, 1, 1, 1}. NURBS basis 
functions has the following properties: (i) non-negativity, (ii) partition of unity, E = 1> ( m ) interpolatory 
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at the end points. As the same function is also used to represent the geometry, the exact representation of 
the geometry is preserved. It should be noted that the continuity of the NURBS functions can be tailored to 
the needs of the problem. The B-spline surfaces are defined by the tensor product of basis functions in two 
parametric dimensions £ and r\ with two knot vectors, one in each dimension as: 

n m 

C (£> V) = E E NiAtWjMVij (35) 

i=l j=l 

where Pjj is the bidirectional control net and Ni tP and Mj q are the B-spline basis functions defined on the knot 
vectors over an m x n net of control points Pjj- The NURBS surface is then defined by: 

C(£,r?) = 7— (36) 

where w(£, rj) is the weighting function. The displacement field within the control mesh is approximated by: 

K> vl,wl, 6%, 6y} = C(£, r]){u oJ , v oJ , w oJ , 9 xJ , 6 yJ }, (37) 
where u o j,v o j,w o j,6 x j,0 y j are the nodal variables and C(£, rj) are the basis functions given by Equation (|36|) , 

4. Shear Locking 

Transverse shear deformations are included in the formulation of Mindlin theory for thick plates. In Mindlin 
theory, the transverse normal to the mid surface of the plate before deformation remain straight but not 
necessarily normal to the mid surface after deformation. This relaxed the continuity requirement on the assumed 
displacement fields. But as the plate becomes very thin, care must be taken in not to violate the following 
relationship 

Vw + 9 = (38) 
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i.e., the shear strain e s must vanish in the domain as the thickness approaches zero. 



Artificial shear correction factor 

Lower order NURBS basis functions, like any other function, suffer from shear locking when applied to thin 
plates. Kikuchi and Ishii 



381 ] introduced an artificial shear correction factor to suppress shear locking in 4-noded 
quadrilateral element. In this paper, we employ the same technique to suppress the shear locking syndrome in 
lower order NURBS basis functions. The modified shear correction factor is given by: 



v P . = v- 



1 + 



fil, 



2n 



1/n 



(39) 



where n, f3 are positive integers, l e is the diameter of the element (=maximum or diagonal length of the element), 
h is thickness of the element and v is the shear correction factor. Here, the shear correction factor obtained 
based on energy equivalence principle as outlined in 



45 



331 ] is used. 



5. Numerical Examples 

In this section, we present the static bending response, free vibration, buckling and flutter analysis of FGM 
plates using a NURBS based finite element method. The effect of various parameters, viz., material gradient 
index n, skewness of the plate ip, the plate aspect ratio a/b, the plate thickness a/h and boundary conditions on 
the global response is numerically studied. The top surface of the plate is ceramic rich and the bottom surface 
of the plate is metal rich. The material properties used for the FGM components are listed in Table [TJ 

Skew boundary transformation 

For skew plates supported on two adjacent edges, the edges of the boundary elements may not be parallel to 
the global axes (x,y,z). In order to specify the boundary conditions on skew edges, it is necessary to use 
the edge displacements (u' , v' , w' Q ), etc., in a local coordinate system (x',y',z') (see Figure ([!])). The element 
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Table 1: Material properties. 



Property 


Aluminum 


Zirconia 


Zirconia 


Alumina 




Al 


Zr0 2 -1 


Zr0 2 -2 


A1 2 3 


E (GPa) 


70 


200 


151 


380 


V 


0.3 


0.3 


0.3 


0.3 


K W/mK 


204 


2.09 


2.09 


10.4 


a/°C 


23 xl0~ 6 


10 xl0~ 6 


10 xlO" 6 


7.2 xlO" 6 


p kg/m 3 


2707 


5700 


3000 


3800 



matrices corresponding to the skew edges are transformed from global axes to local axes on which the boundary 
conditions can be conveniently specified. The relation between the global and local degrees of a particular node 
can be obtained through the following transformation 33 ] 



S = Lgd' 



(40) 



where S and S' are the generalized displacement vector in the global and the local coordinate system, respectively. 
The nodal transformation matrix for a node I on the skew boundary is given by 



cos ip sin ip 
— sin ip cos ip 
1 













cos ip sin ip 



sin ip cos ip 



(41) 



where ip defines the skewness of the plate. 
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5.1. Static Bending 

Let us consider a Al/Zr02 FGM square plate with length-to-thickness a/h = 5, subjected to a uniform load 
with fully simply supported (SSSS) and fully clamped (CCCC) boundary conditions. Four different values 
for the gradient index (n = 0,0.5,1,2) are considered in this study. The plate is modelled with 4, 8, 16, 24 
and 32 control points per side. Tables [2] and [3] summarize the IGA results with quadratic, cubic and quartic 
NURBS elements for SSSS and CCCC boundary conditions. It can be seen that for all polynomial orders, the 
convergence of the results is quite fast. For cubic and quartic NURBS elements, the convergence is almost 
achieved with 16 control points per side. Tab le EH compares the results from the present formulation with other 
approaches available in the literature |46j, |47|, |48|] and a very good agreement can be observed. 
Next, we illustrate the performance of the present isogeometric method for thin plate problems. A simply 
supported and a clamped Al/Zr02-1 square plates subjected to uniform load are considered, while the length- 
to-thickness (a/h) varies from 5 to 10 6 and the gradient index ranges from to 2. Two individual approaches are 
employed: one applied the stabilization technique to eliminate shear locking named S-IGA and the other one, 
normal IGA, without considering any specific technique for shear locking. The plate is modelled using quadratic 
NURBS elements with 13x13 control points. The normalized center deflection w c = 100u> e ^(l-^pa 4 anc ^ the 
normalized axial stress a X x = &xx^z at the top surface of the center of the plate for SSSS and CCCC boundary 
conditions are depicted in Figures [3] - 0J respectively. It is observed that IGA results are subjected to shear 
locking when the plate becomes thin ( a/h > 100). However, the S-IGA results are almost independent of the 



length-to-thickness ratio for thin plates. The same observations also re por ted in 
plates. The results using the S-IGA agree very well with those given in 



49f | for laminated composite 



using the NS-DSG3 element. 



FGM plates under thermo-mechanical loads 

Here, the present isogeometric method is verified on FGM plates subjected to thermo-mechanical loads. A simply 
supported Al/Zr02 plate with aluminium at the bottom surface and zirconia at the top surface is considered. 
The plate with length a = 0.2 m and thickness h=0.01 m is modelled employing quadratic NURBS elements. For 
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Table 2: The normalized center deflection for fully simply supported Al/ZrC>2-l FGM square plate with a/h = 5, subjected to a 
uniformly distributed load p using IGA. *At least 5 control points per side are needed for 1 quartic NURBS element. 



Method 


Number of 




gradient index, n 






Control Points 





0.5 


1 


2 




4 


0.162098 


0.218914 


0.256018 


0.293806 




8 


0.171617 


0.232392 


0.271879 


0.311459 


Quadratic 


16 


0.171649 


0.232439 


0.271935 


0.311520 




24 


0.171651 


0.232441 


0.271938 


0.311523 




32 


0.171651 


0.232442 


0.271938 


0.311523 




4 


0.163329 


0.220598 


0.257991 


0.296053 




8 


0.171658 


0.232452 


0.271950 


0.311536 


Cubic 


16 


0.171651 


0.232441 


0.271938 


0.311522 




24 


0.171651 


0.232442 


0.271938 


0.311523 




32 


0.171651 


0.232442 


0.271938 


0.311523 




5* 


0.172910 


0.234167 


0.273959 


0.313821 




8 


0.171690 


0.232495 


0.272000 


0.311594 


Quartic 


16 


0.171651 


0.232442 


0.271938 


0.311523 




24 


0.171651 


0.232442 


0.271938 


0.311523 




32 


0.171651 


0.232442 


0.271938 


0.311523 
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Table 3: The normalized center deflection for fully clamped Al/ZrC>2-l FGM square plate with a/h = 5, subjected to a uniformly 
distributed load p using IGA. *At least 5 control points per side are needed for 1 quartic NURBS element. 



Method 


Number of 




gradient index, n 






Control Points 





0.5 


1 


2 




4 


0.052510 


0.06818 


0.079306 


0.093421 




8 


0.075831 


0.101036 


0.117946 


0.136557 


Quadratic 


16 


0.076017 


0.101305 


0.118264 


0.136905 




24 


0.076024 


0.101315 


0.118276 


0.136918 




32 


0.076025 


0.101316 


0.118278 


0.136920 




4 


0.056337 


0.07342 


0.085447 


0.100403 




8 


0.076031 


0.101324 


0.118287 


0.136931 


Cubic 


16 


0.076025 


0.101316 


0.118277 


0.136921 




24 


0.076025 


0.101317 


0.118278 


0.136921 




32 


0.076026 


0.101317 


0.118279 


0.136921 




5* 


0.077690 


0.103620 


0.120980 


0.139973 




8 


0.076092 


0.101409 


0.118387 


0.137043 


Quartic 


16 


0.076026 


0.101317 


0.118279 


0.136921 




24 


0.076026 


0.101317 


0.118279 


0.136921 




32 


0.076026 


0.101317 


0.118279 


0.136921 
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Table 4: The normalized center deflection for fully clamped Al/ZrC>2-l FGM square plate with a/h = 5, subjected to a uniformly 
distributed load p using 13 xl3 control points for various boundary conditions. *At least 5 control points per side are needed for 1 
quartic NURBS element. 



Method 


Number of 




gradient index, n 






Control Points 





5 


1 


2 




TO A-OnaHratir 


0.1717 


0.2324 


0.2719 


0.3115 




IGA-Cubic 


0.1717 


0.2324 


0.2719 


0.3115 




IGA-Quartic 


0.1717 


0.2324 


0.2719 


0.3115 


SSSS 


NS-DSG3 [48] 


0.1721 


0.2326 


0.2716 


0.3107 




ES-DSG3 [481 


0.1700 


0.2296 


0.2680 


0.3066 




MITC4 [48J 


0.1715 


0.2317 


0.2704 


0.3093 




A;p-Ritz [47] 


0.1722 


0.2403 


0.2811 


0.3221 




MLPG [46] 


0.1671 


0.2505 


0.2905 


0.3280 




IGA-Quadratic 


0.0760 


0.1013 


0.1183 


0.1369 




IGA-Cubic 


0.0760 


0.1014 


0.1183 


0.1369 




IGA-Quartic 


0.0760 


0.1014 


0.1183 


0.1369 


cccc 


NS-DSG3 [48] 


0.0788 


0.1051 


0.1227 


0.1420 




ES-DSG3 [48] 


0.0761 


0.1013 


0.1183 


0.1370 




MITC4 [48J 


0.0758 


0.1010 


0.1179 


0.1365 




kp-Ritz [47J 


0.0774 


0.1034 


0.1207 


0.1404 




MLPG [46] 


0.0731 


0.1073 


0.1253 


0.1444 
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a/h 

(a) 




(b) 

Figure 3: The normalized center deflection, W c and the normalized axial stress a xx as a function of a/h for a simply supported 
Al/ZrCb-l FGM square plate subjected to a uniform load, p: (a) normalized center deflection and (b) normalized center stress. 
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a/h 



(a) 



IGA, n=0 - + - S-IGA, n=0 -A- IGA, n=0.5 S-IGA, n=0.5 -y- IGA, n=1 SIGA, n=1 — +— IGA, n=2 ■ -X- ■ S-IGA, n=2| 




a/h 



(b) 

Figure 4: The normalized center deflection, w c and the normalized axial stress a xx as a function of a/h for a clamped Al/ZrCb-l 
FGM square plate subjected to a uniform load, p: (a) normalized center deflection and (b) normalized center stress. 
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comparison, a mesh with 13 control points per side is used. At first, we study the behaviour of SSSS Al/Zr02-I 
FGM plate under a uniform mechanical load. Figure © shows the distributions of the normalized axial stress 
through the thickness of the plate computed for different values of the gradient index n. The results are in 
excellent agreement with those given in 



46, 



E c h 3 



481 ] . Figure ([6]) plots the central deflection w c = 100w c 12 (i_^) pa 4 
of the plate with respect to various load parameters, p = pa 4 /(E m h' i ), given in the interval [-14, 0] for different 
values of the gradient index n. It can be seen that the central deflection of the plate linearly increases with 
respect to the load. It is also observed that the central deflection increases with the gradient index k. As 
expected, the metallic plate has the largest deflection, while the ceramic plate has the lowest. Note that the 



results match well with those given in 



5CJ]. Figure (J6j) shows the normalized axial stress at points on the 
vertical line passing through the centroid of simply supported Al/ZrC>2-2 FGM square plate subjected to a 



47|, 



481 ]. It is observed that 



uniform mechanical load p. The results agree well with the results reported in 
for isotropic plates (i.e,n = 0) the axial stress distribution is linear while it is nonlinear for FGM plates. For 
FGM plates, the magnitude of the axial stress at the bottom is less than the one at the top. The maximum 
compressive stress at the top surface of the plate has been obtained for the FGM plate with n = 2, whilst the 
metallic only or ceramic only plates has the minimum tensile stress at the bottom surface. 

FGM plates subjected to thermal loading 

Next, we study the behaviour of FGM plates under the thermal loading. A Al/Zr02-2 FGM plate with length 
a = 0.2m and thickness h = 0.01m is considered. The temperature at the top surface of the plate is varied 
from 0°C to 500° C, while the temperature at the bottom surface is maintained at 20° C. The temperature of 
the stress free state is assumed to be at T a = 0°C. Figure ([7]) depicts the non-dimensional center deflection of 
the plate under the thermal load. This problem was solved by Zhao and Liew, using the element free kp-Kitz 
method 



511 ] and by Xuan et al. 



very well with those given in 



51 



48H u sing NS-DSG3 elements. The results obtained with our IGA method agree 
481 ] . From Figure ([7]), it is observed that the metal plate gives the maximum 



deflection because of its high thermal conductivity. The deflection of the FGM plate with the gradient index 
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Figure 5: The normalized axial stress Un = axxjj^z at points on the vertical line passing through the centroid of a simply supported 
square Al/ZrCb FGM plate subjected to uniform mechanical load p. 



n = 1 is minimum. Generally, it can be seen that the deflections of the FGM plates are much lower than those 
of isotropic plates (i.e.,n = 0), which implies high temperature resistance behaviour of the functionally graded 
plates. 

Now, we investigate the FGM plate under thermo-mechanical loads. The temperature at the top surface of the 
plate is held at 300° C (top surface is assumed to be rich in ceramic) and the temperature at the bottom surface 
(assumed to be rich in metal) is T m =20°C. Figure ([9]) shows the center deflection w c = w c /h of the plate with 
respect to various load parameters p = pa A /(E m h 4 ) given in the interval -14 < p < for different values of 
the gradient index. It can be seen that the central deflection of the plate is completely different from the case 
with a purely mechanical loading (see Figure ©). However, similarly to the case of the FGM plate under pure 
mechanical loading, the center deflection of the plate linearly increases with the load. The metallic phase shows 
the maximum range of deflection changes and the ceramic plate the least. 

Figure ([9]) plots the distribution of axial stress through the thickness of the plate under the uniform mechanical 
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Figure 6: The normalized center deflection w c and the normalized axial stress a xx of a simply supported square Al/ZrCb FGM 
plate subjected to uniform mechanical load p: (a)The normalized center deflection, w c = w c /hwith p — pa /(E m h 4 ) and (b) The 

,2 

normalized axial stress a xx = a xx at points on the vertical line passing through the centroid with load p. 
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Figure 7: The normalized center deflection w c — w c /h of a simply supported square Al/ZrCh FGM plate subjected to thermal load. 



0.3 




Load parameter 



Figure 8: The normalized center deflection w c = w c /h of a simply supported square Al/Zr02 FGM plate subjected to thermo- 
mechanical load. 
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_1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 

Normalized axial stress 

Figure 9: The normalized center deflection w c = w c /h of a simply supported square Al/Zr02 FGM plate subjected to thermo- 
mechanical load. 

load p = -10 6 N/m 2 . Comparing to Figure © for purely mechanical load, it is seen that the maximum 
compressive stress at the top surface of the plate has been obtained for the FGM plate with a gradient index 
n = 1. Again, the metallic or ceramic plate shows the minimum tensile stress at the bottom surface. Note that 



the results match very well with those given in 



51 



48] 



Skew plates 

In this example, we study the behaviour of FGM skew plates under mechanical loads. A simply supported 
Al/Zr02 FGM skew plate with length a = 10m and thickness h = 0.1m is considered. The plate is subjected 
to a uniform mechanical load p = -10 4 N/m 2 . A mesh of quadratic NURBS elements with 17 x 17 control 
points is used for modelling the plate. Figure (|10p shows the distribution of non-dimensional axial stress 
&xx = &xxh 2 / {pa 2 ) through the thickness of the plate for different skew angles with gradient index n = 0.5. It 
can be observed that the axial stresses increase as the skew angle decreases. Similar behaviour can be found for 
gradient indices n = 1 and n = 2 in Figures [11] and [12] respectively. The results obtained by our isogeometric 
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analysis are in a good agreement with those reported in 47(] using the element free kp— Ritz method and the 
results of NS-DSG3 |4g] and ES-DSG3 [52|. 




Non-dimensional axial stress 

,2 

Figure 10: The effect of skew angle on the axial stress a xx = o xx -^j profile for simply supported FGM Al/Zr02 subjected to a 
mechanical load with gradient index n = 0.5. 

5.2. Free flexural vibrations 

In this section, the free flexural vibration characteristics of FGM plates are studied numerically. In all cases, we 
present the non-dimensionalized free flexural frequency defined as, unless otherwise stated: 




where to is the natural frequency, p c , E c are the mass density and Young's modulus of the ceramic phase. Before 
proceeding with a detailed study on the effect of gradient index on the natural frequencies, the formulation 
developed herein is validated against available analytical/numerical solutions in the literature. 
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Non-dimensional axial stress 

Figure 11: The effect of skew angle on the axial stress a xx = o xx -^i profile for simply supported FGM Al/Zr02 subjected to a 
mechanical load with gradient index n = 1.0. 




Non-dimensional axial stress 

Figure 12: The effect of skew angle on the axial stress a xx = cr xx -^ profile for simply supported FGM Al/Zr02 subjected to a 
mechanical load with gradient index n = 2.0. 
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Square plates 

A simply supported AI/AI2O3 FGM square plate with various length-to-thickness ratio is considered. The plate 
is modelled employing quadratic, cubic and quartic NURBS elements with meshes of 8, 14 and 20 control points 
per side. The results obtained from the isogeometric analysis for the first normalized frequency parameter ZJ 
are presented in Tables [5] and [6] for different plate aspect ratios (a/h = 5, 10, 20) and compared with the results 



available in the literature 



0, 



52. 



53 



541 ]. The first normalized frequency (= w '/ 'w analytical) is shown in Figures [T3l 
- Figure (|15p for various a/h ratios (a/h = 5, 10, 20). The results of IGA shown in these figures are computed 
using quadratic NURBS elements with 14 control points per side. It is seen that the results from IGA are 
superior to compared to other methods, irrespective of the material gradient index. 
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Figure 13: The normalized first frequency lj* — ui/uj rB { versus volume fraction exponents for SSSS AI/AI2O3 FGM plate with a/h 
5. The frequency is normalized with the analytical solution. 



Figures [TBI - ()1T|) shows the first four non-dimensionalized frequencies for a Al/Zr02 FGM plate with a/h = 
10 with various boundary conditions (CCCC, SSSS and CFFF). It can be seen that very good agreement is 
obtained with the results available in the literature. For CFFF boundary conditions, the difference between the 
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Figure 14: The normalized first frequency lj* — ui/uj rB { versus volume fraction exponents for SSSS AI/AI2O3 FGM plate with a/h 
10. The frequency is normalized with the analytical solution. 
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Figure 15: The normalized first frequency u* — cj/oj re f versus volume fraction exponents for SSSS AI/AI2O3 FGM plate with a/h 
20. The frequency is normalized with the analytical solution. 
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Figure 16: The non-dimensional frequency oj = ^ versus various volume fraction exponent of Al/ZrCh FGM plate with 

different boundary conditions: (a) Mode 1 and (b) Mode 2. 
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Figure 17: The non-dimensional frequency ui = ^ versus various volume fraction exponent of Al/ZrCh FGM plate with 

different boundary conditions: (a) Mode 3 and (b) Mode 4. 
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IGA and the element free kp— Ritz is ranged from 0.3% to 2.1%, while it is about 2.2 - 3.3% and 1.6 - 2.5% for 
CCCC and SSSS boundary conditions, respectively. Note that all the present results agree very well with those 
given in 48|, [52 ] . 



Table 5: The first normalized frequency parameter ui = uih^Jjf- for fully simply supported AI/AI2O3 FGM square plate for a/h — 
5. 



Method Control points gradient index n 







n 
u 


n ^ 

U.D 


1 


1 

4 


1 n 




Q 
O 


noil oq 


n 1 onc;i 
U.loUOl 


n 1 fiQnn 


n 1 qoco 


n 1 qoqi 


Id A-Oiiadratic 


1 1 

14 


n 01101 


n 1 sn/i k 

U.loU4D 


n 1 fi^n 1 ? 
U.ID0U0 


n 1 ^0^7 


n 1 ^007 




on 


n 01101 


n 1 sn/i /i 

U. 10U44 


n 1 KQn 1 ? 
U.ID0U0 


n 1 qq^7 
u.ioyo 1 


n 1 "?oo7 




S 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 


IGA-Cubic 


14 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 




20 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 




8 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 


IGA-Quartic 


14 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 




20 


0.21121 


0.18044 


0.16303 


0.13957 


0.13227 


ES-DSG3 (20x20) [52] 




0.21218 


0.18114 


0.16351 


0.13992 


0.13272 


DSG3 (16x16) [52] 




0.21335 


0.18216 


0.16444 


0.14069 


0.13343 


MITC4 (16x16) [52] 




0.21182 


0.18082 


0.16323 


0.13968 


0.13251 


HSDT [21] 




0.21210 


0.18190 


0.16400 


0.13830 


0.13060 


kp-Ritz [53] 




0.20550 


0.17570 


0.15870 


0.13560 


0.12840 


Ref. [54] 




0.21120 


0.18050 


0.16310 


0.13970 


0.13240 



Skew plate 

Next, the effects of the skewness of the plate on the free flexural vibration of the FGM plate is studied. A 
Al/Zr02-2 skew FGM plate with length-to-thickness ratio a/h = 10 and various skew angles are considered 
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Table 6: The first normalized frequency parameter ui = uih^Jjf- for fully simply supported AI/AI2O3 FGM square plate for 
10, 20. 



a/h 



Method 




gradient index 


n 









0.5 


1 


4 


10 


IGA-Quadratic (20 points) 


0.05769 


0.04899 


0.4418 


0.03821 


0.03655 


IGA-Cubic (20 points) 


0.05769 


0.04898 


0.04417 


0.03821 


0.03655 


IGA-Quartic (20 points) 


0.05769 


0.04898 


0.04417 


0.03821 


0.03655 


ES-DSG3 (20x20) [52] 


0.05800 


0.04924 


0.04439 


0.03839 


0.03973 


DSG3 (16x16) [52] 


0.05834 


0.04954 


0.04467 


0.03861 


0.03693 


1V111L/4 ^10X10 ) \oZ\ 


U.UO ( ( 


0.049132 


0.04429 


U.U000U 


U.UODDO 


HSDT [21] 


0.05777 


0.04917 


0.04426 


0.03811 


0.03642 


kp-Ritz [53] 


0.05673 


0.04818 


0.04346 


0.3757 


0.03591 


Ref. [54] 


0.05770 


0.04900 


0.04420 


0.03820 


0.03660 


IGA-Quadratic (20 points) 


0.01480 


0.01254 


0.01130 


0.00981 


0.00944 


IGA-Cubic (20 points) 


0.01480 


0.01254 


0.01130 


0.00981 


0.00944 


IGA-Quartic (20 points) 


0.01480 


0.01254 


0.01130 


0.00981 


0.00944 


ES-DSG3 (20x20) [52] 


0.01488 


0.01261 


0.01137 


0.00986 


0.00946 


DSG3 (16x16) [52] 


0.01498 


0.012704 


0.01145 


0.00993 


0.00952 


MITC4 (16x16) [52] 


0.01485 


0.01258 


0.01134 


0.00984 


0.00944 


kp-Ritz [53] 


0.01464 


0.01241 


0.01118 


0.00970 


0.00931 


Ref. [541 


0.01480 


0.01250 


0.01130 


0.00980 


0.00940 
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in this example. The skew plate is modelled using quadratic NURBS elements with 17 x 17 control points. 
The first two non-dimensionalized frequencies are shown in Figures [18] - ()19|) for SSSS and CCCC boundary 
conditions, respectively. For comparison, the results from the element free kp— Ritz |53( ] are also plotted. From 
the figures, it is seen that the IGA gives higher frequencies than the element free kp— Ritz method. It is also 
observed that by increasing the gradient index or decreasing the skew angle, the frequency decreases. In both 
cases, the decrease in the natural frequency can be attributed to the stiffness degradation. In the case of gradient 
index, the stiffness degradation is due to increased metallic volume fraction, while the geometry of the plate is 
a contributing factor in decreasing the frequency when the skew angle decreases. The first eight mode shapes 
of fully clamped Al/Zr02-2 skew plate with skew angle ip = 45° and gradient index n = 0.5 is plotted in Figure 
(HDD. 

Remark 5.1. In the next two sections, the present formulation is extended to study the buckling and flutter 
characteristics of the the FGM plates. For the whole, quadratic NUBRS functions with 17 x 17 control points 
are used, unless otherwise specified. 

5.3. Buckling analysis 

In this section, we present the mechanical and thermal buckling behaviour of functionally graded skew plates. 
The FGM plate considered here consists of aluminum and alumina (see Table Q] for material properties). 

Mechanical Buckling 

The critical buckling parameters are defined for uni- and bi- axial compressive loads as: 



N° b 2 

\ ■* * xxcr 



7T 2 D C 

N*b 2 



A crb = (43) 
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Figure 18: The non-dimensional frequency u = wa 2 ^/ p c /E c for a simply supported Al/ZrC>2-2 skew plate versus various skew angles 
for different values of gradient index with a/h — 10. 
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Mode 1 




(a) mode 1 

Mode 3 




(c) mode 3 

Mode 5 




(e) mode 5 

Mode 7 




(g) mode 7 

Figure 20: First eight mode shapes for a clamped Al/Z: 




(b) mode 2 

Mode 4 




(d) mode 4 

Mode 6 




(f) mode 6 




(h) mode 8 

skew plate with a/h = 10, gradient index n = 0.5 and ip = 45°. 
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where, D c = E c h 3 / '(12(1 — v 2 )). The critical buckling loads evaluated by varying the skew angle of the plate, 
volume fraction index and considering mechanical loads such as uni- and biaxial compressive loads are shown 
in Tables [7] - [S] for two different thickness ratios. The efficacy of the present formulation is demonstrated by 
comparing our results with those in [8J]. It can be seen that increasing the gradient index decreases the critical 
buckling load. It is also observed that the decrease in the critical value is significant for the material gradient 
index n < 2 and that further increase in n yields less reduction in the critical value, irrespective of the skew 
angle. 

Table 7: Critical buckling parameters for a thin simply supported FGM skew plate with a/h — 100 and a/b — 1. 



Skew angle 








Gradient index, 
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1.0830 


Acrb 


2.0002 


1.9999 


0.8980 


0.9017 


0.7660 


0.6303 


0.5415 


15° 


Acru 


4.3946 


4.3946 


1.9716 


1.9716 


1.6752 


1.3800 


1.1868 




-^crb 


2.1154 


2.1154 


0.9517 


0.9517 


0.8086 


0.6652 


0.5716 


30° 


Acru 


5.8966 


5.8966 


2.6496 


2.6496 


2.2515 


1.8607 


1.6032 




-^crb 


2.5365 


2.5365 


1.1519 


1.1519 


0.9788 


0.8044 


0.6905 


45° 


Acru 


10.1031 


10.1031 


4.5445 


4.5445 


3.8625 


3.2234 


2.7964 




Acrb 


3.6399 


3.6399 


1.6863 


1.6863 


1.4330 


1.1774 


1.0103 



Thermal Buckling 

The temperature rise of T m = 5°C in the metal-rich surface of the plate is assumed in the present study. In 
addition to nonlinear temperature distribution across the plate thickness, the linear case is also considered in 
the present analysis by truncating the higher order terms in Equation ([9]) . The plate is of uniform thickness and 
simply supported on all four edges. The critical buckling temperature difference AT cr using two values of the 
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Table 8: Critical buckling parameters for a thick simply supported FGM skew plate with a/h = 10 and a/b = 1. 



Skew angle 








Gradient index, 


n 
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10 






Ref. [8] 


Present 


Ref. [8] 


Present 








0° 


Acru 


3.7374 


3.7307 


1.6892 


1.6793 


1.4198 


1.1632 


0.9999 


A C rb 


1.8686 


1.8654 


0.8449 


0.8397 


0.7099 


0.5816 


0.4999 


15° 


-^cru 


4.0791 


4.0791 


1.8458 


1.8458 


1.5616 


1.2810 


1.1021 




A C rb 


1.9660 


1.9660 


0.8923 


0.8923 


0.7550 


0.6184 


0.5315 


30° 


Acru 


5.3571 


5.3571 


2.4298 


2.4298 


2.0533 


1.6886 


1.4565 




Acrb 


2.3226 


2.3226 


1.0659 


1.0659 


0.9011 


0.7367 


0.6326 


45° 


Acru 


8.5261 


8.5261 


3.8835 


3.8835 


3.2679 


2.7046 


2.3521 




A C rb 


3.1962 


3.1962 


1.5030 


1.5030 


1.2680 


1.0335 


0.8871 



aspect ratio a/b = 1 and 2 with a/h = 10 and for various skew angles is given in Table [9j It can been seen that 
the results from the present formulation are in good agreement with the results available in the literature. The 
decrease in the critical buckling load with the material gradient index n is attributed to the stiffness degradation 
due to the increase in the metallic volume fraction. The thermal stability of the plate increases with the skew 
angle of the plate and the same behavior is observed for other values of gradient index n. It can also be seen 
that the nonlinear temperature variation through the thickness yields higher critical values compared to the 
linear distribution case. 

5.4- Flutter analyses 

In this section, the present formulation is extended to analyse the flutter characteristics of functionally graded 
material plates. Both simply supported and clamped boundary conditions are considered in this study and the 
flow direction is assumed to be at right angles to the plate. Only square plate is considered and the results are 
presented only for a/h = 100. It should be noted that the present formulation is not limited to this alone. In 
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Table 9: Critical buckling temperature AT cr for a thin simply supported FGM skew plate with a/h — 100 and a/b = 1 under linear 
and nonlinear temperature rise through the thickness of the plate. 



a/b 


Skew angle 


Temperature rise 




Gradient index. 


/? 



























Ref. [55] 


Present 


0.5 


1 


5 




0° 


Linear 


24.1951 


24.1912 


9.3787 


5.5207 


3.8987 




Nonlinear 


24.1951 


24.1912 


12.3629 


7.6615 


4.8740 


1 


30° 


Linear 


33.9558 


33.9503 


14.9115 


9.7737 


7.4681 






Nonlinear 


33.9558 


33.9503 


19.6600 


13.558 


9.3399 




60° 


Linear 


123.0974 


123.1172 


65.4519 


48.6271 


40.0647 






Nonlinear 


123.0974 


123.1172 


86.2949 


67.4989 


50.1064 




0° 


Nonlinear 


75.4278 


75.4475 


50.6564 


38.6525 


28.3067 


2 


30° 


Nonlinear 


100.9349 


100.9512 


69.7225 


54.0838 


39.9718 




60° 


Nonlinear 


304.6912 


304.6421 


222.0921 


177.4280 


133.1456 
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all cases, we present the non dimensionalized critical aerodynamic pressure, X cr and critical frequency uj cr as, 
unless specified otherwise: 



_ 2 [Pch 
\L cr — uj cr a y -jj- 

a 3 

X cr = A cr — (44) 

J-Jc 

where D c = 12 ^ c ^a) * s the bending rigidity of the plate, E c , v c are the Young's modulus and Poisson's ratio of 
the ceramic material and p c is the mass density. In order to be consistent with the existing literature, properties 
of the ceramic are used for normalization. 

Table 10: Comparison of critical aerodynamic pressure and coalescence frequency for an isotropic plate with various boundary 
conditions (a/b = 1, a/h = 1 00, v = 0.3,6' = 0). 

Reference Flutter bounds Boundary condition 



1 1 Simply supported Clamped 



Ref. [27 



Present 



A cr 511.11 852.34 

uj cr 1840.29 4274.32 

A cr 511.92 854.88 

uj ct 1844.80 4305.30 



Before proceeding with the detailed study, the formulation developed herein is validated against available results 
pertaining to the critical aerodynamic pressure and critical frequency for an isotropic plate with and without 
a crack. The computed critical aerodynamic pressure and the critical frequency for an isotropic square plate 
with various boundary conditions is given in Table [TUl Next, the influence of boundary conditions on the flutter 
characteristics is studied. For this study, consider a square FGM plate made up of Aluminum- Alumina with 
a/h = 100. Figure (I2ip shows the influence of the boundary conditions on the critical aerodynamic pressure for 
various gradient index. It can be seen that the critical pressure is more for the clamped plate in comparison 
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with that of the simply supported plate as expected. It is also seen that the aerodynamic pressure decreases 
with increase in the gradient index n. However the rate of decrease is high for low values of n. This can be 
attributed to the fact that the stiffness is high for the ceramic plate and minimum for the metallic plate and it 
degrades gradually with increase in the gradient index n. 

900 











ssss 






^^cccc 















4 5 6 
Gradient index, n 



Figure 21: Effect of boundary conditions on the critical aerodynamic pressure A cr 



6. Conclusions 

In this paper, we applied the NURBS based Bubnov-Galerkin iso-geometric finite element method to study the 
static and dynamic response of functionally graded material plates. The first order shear deformation plate 
theory (FSDT) was used to describe the plate kinematics. Of course the present method is not limited to FSDT 
and can easily be extended to higher order plate theories. It is to be noted that with NURBS basis functions, 
geometry could be exactly represented. Although in the present study only simple geometries are considered, 
the only thing that would change is the information pertaining to the geometry represented by the NURBS 
basis functions, when it is applied to model and/or analyze complex geometries. The formulation when applied 
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to thin plates, suffers from shear locking, which is alleviated by employing a modified shear correction factor. 
Numerical experiments have been conducted to bring out the influence of the gradient index, the plate aspect 
ratio and the plate thickness on the global response of functionally graded material plates. From the detailed 
numerical study, it can be concluded that with increasing gradient index n, the static deflection increases, while 
the free flexural vibration, critical buckling load and the flutter frequency decreases. This can be attributed to 
the reduction in stiffness of the material structure due to increase in the metallic volume fraction. 
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